#include <stdio.h>
#include <math.h>
#include <string.h>
#include "slave.h"

__thread_local volatile unsigned long get_reply,put_reply;
__thread_local double A_lhs[BLKY][BLKZ];
__thread_local double A_rhs[3][BLKY+2][BLKZ];
#ifdef _PROFILING
__thread_local volatile unsigned long st,et;
__thread_local volatile unsigned long total_count=0;
__thread_local volatile unsigned long dma_start,dma_end,dma_start_,dma_end_;
extern unsigned long slave_comp_count[64];
extern unsigned long slave_dma_count[64];
#endif
extern double A[2][SIZEX][SIZEY][SIZEZ];
extern int thread_num, iter, N, ntiley,ntilez;
extern int startk[64], startj[64];
extern int reuse_num;
#define alpha_d 0.0876
#define beta_d  0.0765
#ifdef _PROFILING
unsigned long rtc_()
{
  unsigned long rpcc;
  asm volatile("rcsr  %0,4":"=r"(rpcc));
  return rpcc;
}
#endif
void func()
{
   int i,j,k,n,jj,kk,my_id;//N ~ average num tiles that each thread assigned
   int dest,src,my_tiles;
   int cur_blky, cur_blkz;
   my_id = athread_get_id(-1);
      dest = iter % 2 ;
      src = 1-dest;
      n = 0;
      if(my_id<thread_num-1)
         my_tiles=N;
      else
         my_tiles=ntiley*ntilez-N*(thread_num-1);
      //printf("my_tiles:%d\n",my_tiles);
      //my_nums: 
      for(j = startj[my_id]; n < my_tiles && j < SIZEY-2; j += BLKY)
      for(k = startk[my_id]; n < my_tiles && k < SIZEZ-2 ; k += BLKZ)
      {
         int stride,len,bsize;
         ++n;
         //average tile size: <SIZEX, BLXY,BLKZ>
         cur_blkz=BLKZ;
         if(k+BLKZ>SIZEZ-2) //last tile at dim Z
            cur_blkz=(SIZEZ-2)%BLKZ;
         cur_blky=BLKY;
         if(j+BLKY>SIZEY-2) //last tile at dim Y
            cur_blky=(SIZEY-2)%BLKY;
#ifdef _PROFILING
    dma_start=rtc_();
#endif
             //step 1: prelude (start-up)
              get_reply = 0;
              //update low plane 
              stride = (SIZEZ-(cur_blkz))*8;
              len = (cur_blky+2)*(cur_blkz)*8;
              bsize = (cur_blkz)*8;
              athread_get(PE_MODE,&A[src][0][j][k+1]
                 , &A_rhs[0][0][0]
                 , len, (void*)&get_reply
                 , 0, stride, bsize);
             //update middle plane 
             stride = (SIZEZ-(cur_blkz))*8;
             len = (cur_blky+2)*(cur_blkz)*8;
             bsize = (cur_blkz)*8;
             athread_get(PE_MODE,&A[src][1][j][k+1]
                 , &A_rhs[1][0][0]
                 , len, (void*)&get_reply
                 , 0, stride, bsize);
            while(get_reply!=2);
        //step 2: stream in and computation 
         for(i = 1; i < SIZEX-1; i += 1)
         {
           get_reply = 0;
           //update down plane 
           stride = (SIZEZ-(cur_blkz))*8;
           len = (cur_blky+2)*(cur_blkz)*8;
           bsize = (cur_blkz)*8;
           athread_get(PE_MODE,&A[src][i+1][j][k+1]
                 , &A_rhs[(i+1)%3][0][0]
                 , len, (void*)&get_reply
                 , 0, stride, bsize); 
           while(get_reply!=1);
       //step 2.2: computation
#ifdef _PROFILING
       dma_end=rtc_();
       st=rtc_();
#endif              
       for(jj = 1; jj < cur_blky+1; ++jj) 
          for(kk = 1; kk < cur_blkz+1; ++kk)  
             A_lhs[jj-1][kk-1] = alpha_d * (A_rhs[i%3][jj][kk-1]) + \
                        beta_d * (A_rhs[(i-1)%3][jj][kk-1] + A_rhs[(i+1)%3][jj][kk-1] +\
                        A_rhs[i%3][jj-1][kk-1] + A_rhs[i%3][jj+1][kk-1] );
#ifdef _PROFILING
       et=rtc_();
       slave_comp_count[my_id] += ((et-st)); 
       dma_start_=rtc_();
#endif
       //step 3: write back
       put_reply=0;
       stride = (SIZEZ-cur_blkz)*8;
       bsize = cur_blkz*8;
       len = bsize*cur_blky;
       athread_put(PE_MODE,&A_lhs[0][0]
                     , &A[dest][i][j+1][k+1]
                     ,len,(void*)&put_reply,stride,bsize);
       while(put_reply!=1);
       }//of i loop, pipeline loop
    }//of k loop, partitioned loop
#ifdef _PROFILING
    dma_end_=rtc_();
    slave_dma_count[my_id] += (((dma_end_-dma_start_)+(dma_end-dma_start)));
#endif
}
